tau=0.1;    % time step
T=100;
time=tau:tau:T;

type=[];
type(end+1)=-1;      % soma
type(end+(1:2))=0;   % trunk of apical dendrite 
type(end+(1:2))=1;   % C_1
type(end+(1:6))=2;   % C_2
type(end+(1:12))=3;  % C_3
type(end+(1:24))=4;  % C_4
n=length(type);      % the number of compartments

m=max(type)-min(type)+1;  % the number of compartment types

mvr=-60+zeros(1,m);
mvt=[-50,-50+zeros(1,m-1)];
mp=[3,3+zeros(1,m-1)];
mC=[100,100+zeros(1,m-1)];
ma=[0.01,0.01+zeros(1,m-1)];
mb=[5,5+zeros(1,m-1)];
mc=[-55,-35+zeros(1,m-1)];
md=[500,1000+zeros(1,m-1)];
mspike_peak=[50,10+zeros(1,m-1)];

vr=mvr(2+type);
vt=mvt(2+type);
p=mp(2+type);
C=mC(2+type);
a=ma(2+type);
b=mb(2+type);
c=mc(2+type);
d=md(2+type);

spike_peak=mspike_peak(2+type);

ggap=70;
mgapm=ggap*[1 1 1 1 1 1];      % conductance (sensitivity) to j-1 (mother) compartment
mgapd=ggap*[1 1 1 1 1 1];     % conductance (sensitivity) to j+1 (daughter) compartment
gapm=mgapm(2+type);
gapd=mgapd(2+type);


nbrs=zeros(n,3);
nbrs(1,:)=[1 1 2];     % self-gap junctions (no effect), C_0
nbrs(2,:)=[1 2 3];     % C_s, self (no effect), C_0
nbrs(3,:)=[2 4 5];     % C_0, C_1, C_1
nbrs(4,:)=[3 6 7];     % C_0, C_2, C_2
nbrs(5,:)=[3 8 9];     % C_0, C_2, C_2
nbrs(6,:)=[4 12 13];   % C_1, C_3, C_3
nbrs(7,:)=[4 14 15];   % C_1, C_3, C_3
nbrs(8,:)=[5 16 17];   % C_1, C_3, C_3
nbrs(9,:)=[5 18 19];   % C_1, C_3, C_3
nbrs(10,:)=[1 20 21];   % C_s, C_3, C_3
nbrs(11,:)=[1 22 23];   % C_s, C_3, C_3
nbrs(12,:)=[6 24 25];   % C_2, C_4, C_4
nbrs(13,:)=[6 26 27];   % C_2, C_4, C_4
nbrs(14,:)=[7 28 29];   % C_2, C_4, C_4
nbrs(15,:)=[7 30 31];   % C_2, C_4, C_4
nbrs(16,:)=[8 32 33];   % C_2, C_4, C_4
nbrs(17,:)=[8 34 35];   % C_2, C_4, C_4
nbrs(18,:)=[9 36 37];   % C_2, C_4, C_4
nbrs(19,:)=[9 38 39];   % C_2, C_4, C_4
nbrs(20,:)=[10 40 41];   % C_2, C_4, C_4
nbrs(21,:)=[10 42 43];   % C_2, C_4, C_4
nbrs(22,:)=[11 44 45];   % C_2, C_4, C_4
nbrs(23,:)=[11 46 47];   % C_2, C_4, C_4
nbrs(24,:)=[12 24 24];   % C_3, self, self
nbrs(25,:)=[12 25 25];   % C_3, self, self
nbrs(26,:)=[13 26 26];   % C_3, self, self
nbrs(27,:)=[13 27 27];   % C_3, self, self
nbrs(28,:)=[14 28 28];   % C_3, self, self
nbrs(29,:)=[14 29 29];   % C_3, self, self
nbrs(30,:)=[15 30 30];   % C_3, self, self
nbrs(31,:)=[15 31 31];   % C_3, self, self
nbrs(32,:)=[16 32 32];   % C_3, self, self
nbrs(33,:)=[16 33 33];   % C_3, self, self
nbrs(34,:)=[17 34 34];   % C_3, self, self
nbrs(35,:)=[17 35 35];   % C_3, self, self
nbrs(36,:)=[18 36 36];   % C_3, self, self
nbrs(37,:)=[18 37 37];   % C_3, self, self
nbrs(38,:)=[19 38 38];   % C_3, self, self
nbrs(39,:)=[19 39 39];   % C_3, self, self
nbrs(40,:)=[20 40 40];   % C_3, self, self
nbrs(41,:)=[20 41 41];   % C_3, self, self
nbrs(42,:)=[21 42 42];   % C_3, self, self
nbrs(43,:)=[21 43 43];   % C_3, self, self
nbrs(44,:)=[22 44 44];   % C_3, self, self
nbrs(45,:)=[22 45 45];   % C_3, self, self
nbrs(46,:)=[23 46 46];   % C_3, self, self
nbrs(47,:)=[23 47 47];   % C_3, self, self


for swtch=4;
v=zeros(length(time),n);
v(1,:)=vr;
u=0*v;
g_ampa=0+0*v;
I=0+0*v;
st=65/tau;
switch swtch
   case 1
   for i=st:length(time)
       t=time(i)-time(st);
%       g_ampa(i,3)=100*t*exp(-2*t);
   end;
   case 2
   for i=st:length(time)
       t=time(i)-time(st);
       g_ampa(i,1)=300*t*exp(-2*t);
   end;
   case 3
   for i=st:length(time)
       t=time(i)-time(st);
       g_ampa(i,1)=300*t*exp(-2*t);
   end;
   I(:,:)=70; % Input current starts high for all compartments, so action potential propagates to distal dendrites.

   case 4
   for i=st:length(time)
       t=time(i)-time(st);
       g_ampa(i,24)=300/2*t*exp(-2*t);
       g_ampa(i,25)=300/2*t*exp(-2*t);
   end;
   case 5
   for i=st:length(time)
       t=time(i)-time(st);
       g_ampa(i,24)=300*t*exp(-2*t);
   end;
   case 6
   for i=st:length(time)
       t=time(i)-time(st);
       g_ampa(i,24)=300*t*exp(-2*t);
   end;
   I(:,:)=60;
   otherwise
end;

for i=1:length(time)-1    
 dv=p.*(v(i,:)-vr).*(v(i,:)-vt) -u(i,:) + I(i,:) + g_ampa(i,:).*(0-v(i,:)) + ...
     gapm.*(v(i,nbrs(:,1))-v(i,:)) + ...
     gapd.*(v(i,nbrs(:,2))+v(i,nbrs(:,3))-2*v(i,:));
 v(i+1,:)=v(i,:)+tau*dv./C;
 u(i+1,:)=u(i,:)+tau*(a.*(b.*(v(i,:)-vr)-u(i,:)));
 fired =find(v(i+1,:)>spike_peak);
 v(i,fired)=spike_peak(fired);
 v(i+1,fired)=c(fired);
 u(i+1,fired)=u(i+1,fired)+d(fired);
end;

switch swtch
   case 1
       %subplot(4,3,1)
       plot(time,v(:,1),time,v(:,3),time,v(1,1)+4*(v(:,1)-v(1,1)),time,-2+v(1,1)+0.05*g_ampa(:,3));
       axis([50 100 -63 -53])
   case 2
       subplot(1,5,1)
       ind=[1 2 3 4 6 12 24 25];
       for i=1:length(ind)
           ii=dilute(v(:,ind(i)),-40,10);
           plot(time(ii),(i-1)*(80)+v(ii,ind(i))); hold on
       end;
       ii=dilute(g_ampa(:,1),40,10);
       plot(time(ii),-99+0.5*g_ampa(ii,1),'r')
       hold off;
       xlim([50 100])
       ylim([-100 600]);
   case 3
       subplot(1,5,2)
       ind=[1 2 3 4 6 12 24 25];
       for i=1:length(ind)
           ii=dilute(v(:,ind(i)),-40,10);
           plot(time(ii),(i-1)*(80)+v(ii,ind(i))); hold on
       end;
       ii=dilute(g_ampa(:,1),40,10);
       plot(time(ii),-99+0.5*g_ampa(ii,1),'r')
       hold off;
       xlim([50 100])
       ylim([-100 600]);
       axis off;
   case 4
       subplot(1,5,3)
       ind=[1 2 3 4 6 12 24 25];
       for i=1:length(ind)
           ii=dilute(v(:,ind(i)),-40,10);
           plot(time(ii),(i-1)*(80)+v(ii,ind(i))); hold on
       end;
       ii=dilute(g_ampa(:,24),40,10);
       plot(time(ii),-99+6*(80)+0.5*g_ampa(ii,24),'r')
       ii=dilute(g_ampa(:,25),40,10);
       plot(time(ii),-99+7*(80)+0.5*g_ampa(ii,25),'r')
       hold off;
       xlim([50 100])
       ylim([-100 600]);
       %axis off;
   case 5
       subplot(1,5,4)
       ind=[1 2 3 4 6 12 24 25];
       for i=1:length(ind)
           ii=dilute(v(:,ind(i)),-40,10);
           plot(time(ii),(i-1)*(80)+v(ii,ind(i))); hold on
       end;
       ii=dilute(g_ampa(:,24),40,10);
       plot(time(ii),-99+6*(80)+0.5*g_ampa(ii,24),'r')
       hold off;
       xlim([50 100])
       ylim([-100 600]);
       %axis off;
   case 6
       subplot(1,5,5)
       ind=[1 2 3 4 6 12 24 25];
       for i=1:length(ind)
           ii=dilute(v(:,ind(i)),-40,10);
           plot(time(ii),(i-1)*(80)+v(ii,ind(i))); hold on
       end;
       ii=dilute(g_ampa(:,24),40,10);
       plot(time(ii),-99+6*(80)+0.5*g_ampa(ii,24),'r')
       hold off;
       xlim([50 100])
       ylim([-100 600]);
       axis off;
   otherwise
end;

end;

